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Abstract 

Significant differences exist among literature for thermal conductivity of various systems com- 
puted using molecular dynamics simulation. In some cases, unphysical results, for example, neg- 
ative thermal conductivity, have been found. Using GaN as an example case and the direct non- 
equilibrium method, extensive molecular dynamics simulations and Monte Carlo analysis of the 
results have been carried out to quantify the uncertainty level of the molecular dynamics methods 
and to identify the conditions that can yield sufficiently accurate calculations of thermal conduc- 
tivity. We found that the errors of the calculations are mainly due to the statistical thermal 
fluctuations. Extrapolating results to the limit of an infinite-size system tend to magnify the er- 
rors and occasionally lead to unphysical results. The error in bulk estimates can be reduced by 
performing longer time averages using properly selected systems over a range of sample lengths. If 
the errors in the conductivity estimates associated with each of the sample lengths are kept below 
a certain threshold, the likelihood of obtaining unphysical bulk values becomes insignificant. Using 
a Monte-Carlo approach developed here, we have determined the probability distributions for the 
bulk thermal conductivities obtained using the direct method. We also have observed a nonlinear 
effect that can become a source of significant errors. For the extremely accurate results presented 
here, we predict a [0001] GaN thermal conductivity of 185 W/K • m at 300 K, 102 W/K • m at 
500 K, and 74 W/K • m at 800 K. Using the insights obtained in the work, we have achieved a 
corresponding error level (standard deviation) for the bulk (infinite sample length) GaN thermal 
conductivity of less than 10 W/K • m, 5 W/K • m, and 15 W/K • m at 300 K, 500 K, and 800 K 
respectively. 
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I. INTRODUCTION 



The thermal transport property of semiconductor nanostructures is becoming an increas- 
ingly prominent focus of research [THS]. In electronics applications, the decrease in feature 
sizes to the nanometer scale has resulted in significant increases in heat generation. This 
trend has placed more stringent demands on the ability of devices to dissipate heat efficiently 
to the surrounding environment [1J . By contrast, thermoelectrics [6] require a low thermal 
conductivity for efficient performance. In materials engineered with nanometer-scale fea- 
tures, thermal conduction can be greatly controlled by density of interfaces and defects. 
At this "phonon engineering" level, experimental measurements are quite challenging, and 
the path to improved properties is not always clear. For example, it is often not possible 
in experiments to distinguish the contributions of individual defects to thermal resistance. 
To compliment experimental studies, atomic-scale simulation is beginning to play a critical 
role in achieving greater fundamental understanding and identifying improved strategies for 
tailored thermal properties. 

The two most common approaches for computing thermal conductivity based on 
molecular-dynamics (MD) simulation are the Green-Kubo method [7HT3] and the "direct 
method" [HH23]. The Green-Kubo approach involves an equilibrium MD simulation, and 
the thermal conductivity is determined from the time-dependence of the current-current 
correlations functions. In the direct method, a heat current J is applied and the time- 
averaged temperature gradient dT/dx is computed. The thermal conductivity k is then 
obtained from Fourier's law 



Both the direct and Green-Kubo approaches require long simulations (e.g. at least 1 ns) 
to reduce the uncertainty due to thermal fluctuations. For the direct method, another 
difficulty encountered is that the computed thermal conductivity depends strongly on the 
system length L along the propagation direction, which is typically limited to at most a few 
hundred nanometers. This means that for perfect bulk crystals the phonon mean free path is 
comparable to the system size and transport occurs in a partially ballistic regime [T71 \18\ I2H- 
[23] . It follows from kinetic theory that the computed values of k are smaller than that of a 
true bulk system. To obtain values that can be meaningfully compared with experiments, 
it is therefore necessary to perform several simulations for different cell lengths, and then 
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extrapolate to the infinite-size limit. Simple theoretical considerations [ToTfTTt I2T] based on 
the assumption that scattering due to the finite-size simulation cell acts independently from 
other scattering mechanisms suggests that the computed value k depends on the system 
length L as: 

I I a 

- = — + - (2) 

K /too *-> 

where /too is the extrapolated thermal conductivity, and a is a length-independent coeffi- 
cient. Eq. [5] can be recognized as a Matthiessen rule |24] . Some previous MD simulations [21] 
appeared to have shown agreement with the linear dependence of 1/k on 1/L assumed in 
Eq. [2} However, there are some instances where the use of Eq. [2] has failed to give reasonable 
results for /too- For example, a recent use of Eq. [2] resulted in an unphysical prediction for 
/too of silica [23 J. It is unclear whether this result is due to a failure of Eq. [2] to describe the 
length dependence of k, or to inaccuracy caused by thermal fluctuations. Previous direct 
method calculations [T4TI23] have typically been performed for small system dimensions (e.g. 
from hundreds to thousands of atoms [TU [16j fl8H20| [23] ) and short simulation times (e.g. 
from a few hundred ps up to a few ns [T4T - [T6"l IT8H22] ). Studies have shown that these small, 
short-time simulations can result in statistical errors of ±10% or more [2"T| |2"2"] . Furthermore, 
extrapolating to the infinite-size limit using Eq. [2] tends to magnify statistical errors and 
/too is very difficult to determine accurately. It should also be recognized that the errors 
associated with the Green-Kubo method are often in the 15% — 20% range or above, de- 
pending on the total simulation time[HHTT]. As a result, the utility of atomistic simulations 
in studying defects (e.g., vacancy, interstitial, etc.) whose effects on thermal conductivity 
are below 20% has been limited. 

In this work, we apply large computer clusters to perform direct method calculations of 
the thermal conductivity of perfect GaN wurtzite crystals along the [0001] direction. GaN 
is of interest due to its desirable optoelectronic properties and its ability to integrate with 
existing silicon structures. Notable applications, such as laser diodes and high electron 
mobility transistors [2"BTI3T)] . operate at high current and power densities and thus accurate 
estimation of heat dissipation due to conduction is crucial. The simulations use expanded 
sample space, relatively large systems (up to 60,000 atoms) and extremely long times (up 
to 40 ns). This reduces the statistical error of each sample at a specific length to the range 
of ±0.5 — 2.0%. To quantify the error in the extrapolated value /too, we have developed 
and applied a Monte-Carlo algorithm to compute statistical distributions of /too- We find 

4 



that is very sensitive to numerical uncertainty in the data points used to fit Eq. [2} 
In particular, we find evidence that the occasionally discovered unphysical results origin 
from numerical uncertainty in each simulated thermal conductivity value. While numerical 
uncertainties are clearly an important issue in direct-method calculations, the high quality 
of the data reported here also enables a direct investigation of the ability of Eq. [2] to describe 
the length-dependence of simulated data. In contrast to many other studies, we find strong 
evidence that some nonlinear effects are present which cannot be described by Eq. [2j Finally, 
we also have assessed the dependence of the simulated thermal conductivity on various 
parameters including cross-sectional area, size of the hot and cold reservoirs, the heat flux, 
and the thermal expansion freedom. We compare the results with those obtained using 
selected Green-Kubo calculations and with previously published data from both atomic- 
scale simulations and experiments. 

II. METHOD 

A. Interatomic potential 

There are at least two different Tersoff potentials [311 E2] and two different Stillinger- Weber 
(SW) potentials [33H35] with parameters developed for GaN. Some thermal conductivity 
studies have been reported from simulations using SW potentials, including transport in 
bulk[71 [36] and nanowires[3~T]. Here, we use the SW potential with parameters developed 
by Bere and Serra[3U [35]. We will compare the results found using this SW potential with 
previously published results [7J [36] . 

To evaluate the suitability of the SW potential [3^1 133] for thermal transport calculations, 
we first apply it to compute the dynamical properties including the dispersion relations, 
vibrational density of states (DOS), and heat capacity using the established theory [38]. The 
results are compared with corresponding experiments |39lHT] in Figs. [T] - [3j 

Fig. [I] shows phonon dispersion relations along the [0001] direction. It can be seen that 
the SW results of the longitudinal acoustic (LA) branch is in very good agreement with 
those obtained from the Raman scattering and inelastic X-ray scattering (IXS) experiments 
[39] . By contrast, the SW potential significantly underestimates the longitudinal optical 
(LO) branch. 
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FIG. 1: Comparison between SW calculations and experimental phonon dispersion relations for 
bulk GaN along the [0001] direction. 

Fig. [2] shows the vibrational DOS data. Consistent with the dispersion curves in Fig. [T], 
the calculated vibrational DOS in Fig.|2]is seen to be in reasonable agreement with the time- 
of-flight neutron spectroscopy experiments [40J for the lower- frequency modes which include 
the acoustic modes. By contrast, there is a substantial difference between the computed 
and experimental frequencies for the high-frequency optical modes. It is noted that the SW 
potential does not include the long-range electrostatic interaction, which is known to be 
responsible for splitting the longitudinal optical (LO) and transverse optical (TO) phonon 
branches in polar materials. However, we also note that the interatomic potentials with the 
electrostatic interactions [42J may also fail to describe both acoustic and optical phonons 
at the same time. 

Fig. [3] shows specific heat C p . It can be seen that the SW predictions are in reasonable 
agreement with the experimental data[H] at low temperatures. At higher temperatures, 
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FIG. 2: Comparison between SW calculations and experimental data of DOS for bulk GaN. 

the calculations tend to underestimate the experimental values. This is consistent with a 
significant overestimation of the optical phonon frequencies by the SW potential. 

In summary, the dynamical properties of the SW potential for GaN exhibits reasonable 
behavior when compared to experiment. For low-frequency modes the agreement is much 
better than for high-frequency optical modes. 

B. Computational cell 

The equilibrium GaN has a wurtzite hexagonal crystal structure. The experimental lattice 
constants are a = 3.19A, c = 5.19A, and an internal displacement u = 0.377|l3]. With the 
SW interatomic potential used here, the zero temperature lattice constants are a = 3.19A, 
c = 5.20A, and u = 0.375. The computational supercell is aligned so that the x—, y—, 
and z— coordinates correspond respectively to [0001], [1100], and [1120] directions. Of the 
hexagonal structure, we can define a unit orthogonal cell whose dimensions in the x—, y—, 
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FIG. 3: Comparison between SW calculations and experimental data[41j for specific heat of bulk 
GaN. 

and z— directions are respectively c, 2-a-cos (7r/6), and a. Along the [0001] x- direction of 
heat propagation, the number of unit cells rti is chosen from the range 150 < n\ < 500, 
which approximately corresponds to the range between 770 A and 2550 A. The numbers 
of unit cells in the [1100] y- and [1120] z- directions, n<i and 77,3, are chosen in the range 
2 x 3 < 112 x n-z < 6 x 10 so that the corresponding cross-sectional area ranges between 
approximately 106 A 2 and 1060 A 2 . 

Initial crystals were created by assigning atom positions according to prescribed crystal 
lattice. Two types of initial crystals were used to explore the effect of thermal expansion. 
The first one was assigned the known lattice constants at the zero temperature. The other 
one was assigned the thermally-expanded lattice constants at the simulated temperature. 
To determine the thermally-expanded lattice parameters, a simulation in the NPT (constant 
atom number, pressure, temperature) ensemble was performed for a bulk GaN crystal for 
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a total of 3 x 10 4 MD steps with a time step size At = lfs. The cell dimensions were 
taken from a time-average over the final 10 4 MD steps. The averaged cell sizes gave well- 
converged thermal-expansion strains (with respect to the equilibrium sizes at K) of about 
0.00194, 0.00287, and 0.00439 at 300 K, 500 K and 800 K temperatures respectively. Once 
initialized according to the designated dimension, the volume of the crystal was conserved 
during the subsequent NVE (constant atom number, volume, and energy) thermal transport 
simulations. 

The initial average temperature was established by assigning velocities to atoms according 
to Boltzmann distribution and zero total linear momentum[20] HH US]. Considering that 
half of the initial kinetic energy is transferred to potential energy due to the equipartition 
of energy and the system energy does not change under the NVE condition, we therefore 
assigned the initial velocities consistent with double the desired temperature. The thermal 
transport simulation is started immediately without the conventional long NPT or NVT 
simulation to establish the initial temperature. An advantage of this approach is that once 
steady-state is reached, the average temperature of the system matches exactly the desired 
temperature, thereby eliminating one possible source of errors due to uncertainty of the 
initialized system energy typically seen in the NPT or NVT runs. 

C. Heat conduction algorithm 

Fig. [4] shows a schematic illustration of the computational cell using an x-y projection. 
As can be seen from Fig. |4j the GaN crystal is composed of a stack of Ga and N bi-plane 
units in the x- direction, with a small separation distance between the two (Ga and N) 
planes within each unit (highlighted with a white area), and a wide separation distance 
between units. To avoid effects of free surfaces and to simulate bulk configurations, our 
model employs periodic boundary conditions in all three coordinate directions. The "direct 
method" requires the creation of a "hot" and a "cold" region. In Fig. |4j the hot region is at 
the far left of the cell, and the cold region is near the middle of the cell. With appropriate 
choices of system length and the boundary positions and width of the hot and cold regions, 
we can ensure that the hot and cold regions are exactly identical, and that the left side of 
the cold (or hot) region is exactly symmetric to its right side up to another hot (or cold) 
region in the periodic image. 



9 



black and white circles distinguish Ga and N 




► X 

FIG. 4: Schematic illustration of the x-y projection of the computational cell. Boundaries of 
periodic cell and heat / cold regions are carefully set at the middle between widely separated planes 
as highlighted respectively by dashed frame and dark shaded areas. 

The heat flux in the direct method can be generated using either constant temperature [i~4T 
\TE\ |36| [37] or constant flux [i~9] - |23"] control for the hot and cold regions. Using constant 
temperature control methods such as Nose-Hoover algorithm[46j and velocity rescaling, we 
found that the thermal conductivity data can vary as a function of the parameters, e.g. 
Nose-Hoover mass or rescale frequency, that affect the "equilibration time" the algorithm 
used to control the temperature, and the heat added to the hot region or removed from the 
cold region is not continuous (jumps from positive to negative values between time steps). 
Here we used the approach of Ikeshoji and Hafskold|44J to create a constant heat flux. In 
constant heat flux simulations, a constant amount of energy is added to the hot region and 
exactly the same amount of energy is removed from the cold region (while preserving linear 
momentum) at each MD time step using velocity rescaling. In our simulations where SW 
potential and an MD time step of 1 fs were used, the total energy was conserved extremely 
well (e.g., the energy drift was about 8 x 10~ 12 eV/ns ■ atom at 300 K). 

Simulations were performed at temperatures of 300 K, 500 K and 800 K. To generate 
extremely accurate results, the duration of simulations was chosen to be at least 24 ns at 
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300 K and at least 44 ns at 500 K and 800 K. To compute the temperature profile, 100 cells 
were created along the x— direction. A temperature averaged over a designated number of 
time steps was calculated for each of the cells. The temperature profile and the input heat 
flux were used to calculate the thermal conductivity using Fourier's Law in Eq. [I] 

III. RESULTS 

A. Statistics and uncertainty of direct method thermal conductivity calculations 

To ensure that the system is at steady state before computing the time-averaged tem- 
perature profile, we neglect the first 4 ns of simulation time after the heat source and sink 
are switched on. For a system with n\ = 500, n 2 = 3, and n 3 = 5 (a total of 60000 atoms) 
at thermally-expanded dimensions, a simulation was carried out using a heat flux of 0.0015 
eV/ps ■ A 2 , an average temperature of 300 K, and a heat source width of 60 A. We show 
in Fig. [5] two temperatures computed for cells 230 A left and right of the cold source. The 
data in Fig. [5^ a) was averaged over 4 ps intervals, whereas the data in Fig. ^b) is the 
running- average over all time steps after the initial 4 ns. The fact that there is no apparent 
drift in the temperature and the two mean temperatures measured at equal distances from 
the cold source are essentially identical in Fig. [5]^a), indicates that the system has attained 
steady state during this period. However, Fig. [5^a) indicates a significant scatter of data. 
The running average in Fig. ^b) shows that if the temperature was averaged over a large 
enough number of time steps, this scatter can be reduced to a negligible value. Note that 
the minimum number of time steps needed to get accurate results was found to increase 
when temperature was increased or when materials with lower thermal conductivity were 
simulated. 

To elucidate the dependence of the computed temperature profile on the total averaging 
time, we show in Fig. [6] temperature profile data points computed for averaging times of 0.5 
ns, 1.0 ns, 4.0 ns, and 20.0 ns. For a heat flux of 0.0015 eV/ps ■ A 2 , a linear region always 
exists in the temperature profile away from the heat source and sink. The temperature 
gradient was determined by fitting to only the linear regions which were taken to occupy the 
middle half of the length between the heat source and sink. The fitted linear functions are 
shown in Fig. [6] using the thick gray lines along with the computed temperature gradients. 
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FIG. 5: Temperature convergence curves. 

For averaging times less than 1 ns, the profiles show significant thermal fluctuations. In 
addition, the temperature gradients measured at both sides of the cold region did not match 
exactly, indicating statistical errors. By contrast, for averaging times of 4 ns and greater, 
the temperature profiles became progressively smoother. In addition, the difference between 
the two temperature gradients became very small when the averaging time reached 20 ns. 
The evolution of the temperature gradient as a function of averaging time can be more 
clearly seen from the running averages of the two temperature gradients shown in Fig. [7} 
It can be seen that between 5 and 20 ns time, the left and the right temperature gradients 
were considerably different and there is a decreasing trend with time in the difference of the 
temperature gradients. Note that our system was initialized with a uniform temperature 
distribution, and the nominal temperature gradient was uniformly zero at the start of the 
simulation. The difference in the left and the right temperature gradients, however, occurred 
due to system transition from the initial state that includes some randomness of the initially 
assigned atom velocities. Fig. [7] again indicates that for this particular simulated case, an 
averaging time of greater than 20 ns is required to reduce thermal fluctuation to a negligible 
level. Note that much longer averaging time would be needed if the simulation was at 800 
K. 
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FIG. 6: Time evolution of temperature profiles. 
B. Thermal conductivity calculation errors 

Using the same system and the same heat source width, we further analyze the numerical 
errors of the calculated thermal conductivity data k. After discarding the first 4 ns of 
simulation time, the remaining simulation time was uniformly divided into TV = 20 segments. 
The thermal conductivity k was computed from the temperature profiles averaged over each 
of the time segments as well as the running temperature profiles averaged over all time 
steps. The results obtained at 300 K temperature, 0.0015 eV/ps ■ A 2 heat flux and those 
at 800 K temperature, 0.0008 eV/ps ■ A 2 heat flux are shown respectively in Figs. |8|a) and 
fflh), where the data points and lines represent respectively the short-time average and the 
running average values. 

It can be seen from Fig. [8] that the short-time-averaged thermal conductivity data is 
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FIG. 7: Running average of the absolute value of temperature gradients measured on the left and 
the right side of the cold region. 

very scattered. At 300 K, k ranges from about 85 W/K ■ m to 100 W/K ■ m. At 800 K, 
k ranges from under 45 W/K ■ m to about 70 W/K • m. Notice that the scatter is greater 
for simulations at higher temperatures especially considering that the time segments at 300 
K are of duration 1 ns whereas those at 800 K are of 2 ns. The running averages indicate, 
however, that the error is sharply reduced by increased averaging times. 

The best estimate of the thermal conductivity, Ki, for a system of length of Lj, is deter- 
mined from the temperature profiles averaged over the entire time of simulation excluding 
the initial 4 ns. The standard deviation of the calculated thermal resistivity (inverse of 
thermal conductivity) can be calculated from the short-time averaged thermal conductivi- 
ties. Here resistivity rather than conductivity is considered since it is linearly related the 
term dT/dx, which is the error source in Eq. [I] (i.e. given that in the Ikeshoji and Hafskold 
algorithm the heat current J is exact to machine precision) . We first compute M short-time 
averaged thermal conductivities k^i, k^, k^a/"- The standard deviation of the resistivity 
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FIG. 8: Statistics of thermal conductivity calculations, 
of the short-average data, fjj s , is expressed as, 



where 



1 



spN f J. 1_ 



(3) 



(4) 



is the short-time averaged thermal resistivity. Because dT/dx is in the numerator, the long- 
time averaged thermal resistivity is an average of the short-time thermal resistivity. The 
standard deviation of the long-time averaged thermal resistivity 1//%, therefore, is given by 



(5) 



To provide guidance for the choice of averaging time in thermal conductivity calculations, 
we show in Fig. [9] the standard deviation of thermal resistivity as a function of averaging 
time using the thermally-expanded 500 x 3 x 5 sample simulated at 800 K. It can be seen 
that the standard deviation was reduced rapidly when the averaging time was increased 
from 4 ns (total simulation time 8 ns) to 20 ns (total time 24 ns), but the rate of decrease 
becomes small when the averaging time was further increased. From Eq. [5j it is apparent 
that the standard deviation rjj should depend on J\f as J\f~z, which is consistent with the 
results in Fig. [9j Based on an estimate equation, 



AKi « K; ■ Oi 



(6) 



15 



0.0030 



t*H 0.0025 



■i 0.0020 

09 

• 7-H 

!*) 

<u 
u 

% 0.0015 

e 

o 

■ — 

"§ 0.0010 
> 

"2 0.0005 

-a 

c 

» 0.0000 



- 1 r 

o 



T 1 1 1 - 



800 K, 500x3x5 sample 



o o 



o o 



o o 



-I l_ 



4 8 12 16 20 24 28 32 36 40 

average time (ns) 



FIG. 9: Standard deviation of thermal conductivity as a function of averaging time for a selected 
example of simulation. 

the standard deviation of resistivity (<7j) needs to be less than 0.0003 K ■ m/W for the 
standard deviation of conductivity (A«j) to be less than 1.0 K • m/W for the 800 K case. 
This means that the averaging time needed for the calculation is significantly longer than 
50 ns. 



C. Effects of cross-sectional area, source width, and heat flux 

We have also explored the dependence of the computed thermal conductivity on cross- 
sectional area, source width, and the magnitude of the heat flux. The dependence of the 
computed thermal conductivity on these parameters is shown respectively in Figs. |To|a) 



- fLOTc). In Fig. 10 a), it is apparent that the dependence on cross-sectional area is quite 



weak, although the conductivity is slightly higher for the smallest cross-sectional area of 2 x 3. 



Fig. 10 b) shows that the computed thermal conductivity is independent of the magnitude 



of the heat flux, indicating linear behavior within this range. In addition, we see from 
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FIG. 10: Dependence of calculated thermal conductivity upon cross-sectional area, heat flux, and 
source width. 



Fig. |10[c) that the results do not depend on the width of the hot and cold sources. These 
observations are in agreement with previous results pi], but are established more definitively 
here due to the extremely small statistical error. We also found that for a given averaging 
time, simulations for larger cross-sectional area obviously tend to produce smaller standard 
deviations due to the larger number of atoms that is averaged in each cell. It is therefore the 
case that choosing a small cross-sectional area does not decrease the computational costs as 
dramatically as might be expected, since small cross-sectional areas require longer averaging 
times. 

In summary, we have identified a wide range of values for the cross-sectional area, source 
width, and heat flux that do not affect the computed thermal conductivity. However, we 
will see in the next section that there do appear to be some instances at higher temperatures 
and for longer sample sizes where results can depend on the magnitude of the heat current 
J. This behavior, as we will show in the next subsection, appears to result in deviations 
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TABLE I: Thermal conductivity m, standard deviation of thermal resistivity Oj, and output tem- 
perature Ti obtained at different sample length L{ (i = 1, 2, 8) but fixed input temperature of 
300 K, heat flux of 0.0015 (eV/ps - A 2 ), and averaging time of 20.0 ns including thermal expansions. 
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from the length dependence predicted by Eq. [2] 
D. GaN thermal conductivity 

Having rigorously established the inherent statistical error as a function of simulation 
time, we have computed GaN thermal conductivities at temperatures of 300 K, 500 K, and 
800 K. Simulations were performed for 8 different sample sizes at each temperature, covering 
a wide range of sample length (L, = n\ from 150 to 500 cells), to explore extrapolation using 
Eq. [2] Each system had a cross-sectional area of 3 x 5 unit cells, and a source width of 60 A. 
At each simulated temperature, thermal expansion of the supercell was included. However, 
to explore whether thermal expansion plays an important role, we also computed the thermal 
conductivity at 300 K and 500 K with the lattice parameters fixed at their K values. In 
Tables [I] -[V] we present the results for thermal conductivity Ki, standard deviation of thermal 
resistivity <7j, and average temperature Ti for each of the system length Li (i=l,2,...,8). Other 
parameters used in these simulations are given in the table captions. 

Tables [T]- [V] show that the average temperature, which is the point at which the tempera- 
ture gradient (and hence the thermal conductivity) was calculated, is very consistent among 
different runs (within 0.5 K from the mean output temperature for each series). Based upon 
the approximate relation, Eq. [6j Tables [l| - |V| also show that the averaging time used (20.0 
ns for 300 K series and 40.0 ns for 500 K and 800 K series) results in standard deviations 
of thermal conductivity e» of less than 1.0 W/K ■ m at 300 K and 500 K and less than 
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TABLE II: Thermal conductivity m, standard deviation of thermal resistivity cii, and output 
temperature Tj obtained at different sample length Lj (i = 1, 2, 8) but fixed input temperature of 
500 K, heat flux of 0.0012 (eV/ps-A 2 ), and averaging time of 40.0 ns including thermal expansions. 



i 


1 


2 3 


4 


5 


6 


7 


8 


Li(cells) 


150 


200 250 


300 


350 


400 


450 


500 


«i(W/K-m) 


29.36 


34.05 38.77 


43.72 


47.70 


52.31 


55.97 


60.42 


<Ji{K-m/W) 


0.000654 


0.000463 0.000376 


0.000275 


0.000383 


0.000197 


0.000184 


0.000199 


Ti 


501.00 


500.68 501.02 


501.09 


501.23 


501.34 


501.14 


501.24 






= 101.64 W/K • m, 


Koo,MC = 


101.80 ± 3. 


88 W/K • 


m 





TABLE III: Thermal conductivity Kj, standard deviation of thermal resistivity <jj, and output 
temperature Tj obtained at different sample length Lj (i = 1, 2, 8) but fixed input temperature of 
800 K, heat flux of 0.0008 (eV/ps-A 2 ), and averaging time of 40.0 ns including thermal expansions. 



Li(cells) 


150 


200 250 300 350 400 


450 


500 


Ki(W/K-m) 


20.37 


24.85 28.15 35.40 39.94 45.72 


51.96 


55.26 


<Ji{K-m/W) 


0.001265 


0.001235 0.000861 0.000931 0.000606 0.000668 


0.000511 


0.000440 


Ti 


801.21 


801.38 801.35 801.40 801.64 801.59 


801.86 


801.80 




Koo,ext 


= 71.88 W/K • m, k^ mc = 73.76 ± 12.73 W/K • 


m 





TABLE IV: Thermal conductivity Kj, standard deviation of thermal resistivity Oj, and output 
temperature Tj obtained at different sample length Lj (i = 1, 2, 8) but fixed input temperature 
of 300 K, heat flux of 0.0015 (eV/ps-A 2 ), and averaging time of 20.0 ns without thermal expansions. 



i 


1 


2 3 


4 


5 6 


7 


8 


Li(cells) 


150 


200 250 


300 


350 400 


450 


500 


Ki(W/K-m) 


42.39 


51.35 58.92 


66.16 


74.83 79.15 


85.70 


90.43 


(Ti(K-m/W) 


0.000432 


0.000280 0.000239 


0.000193 


0.000129 0.000132 


0.000114 


0.000113 


Ti 


300.59 


300.74 300.76 


300.88 


300.77 300.93 


300.77 


301.04 




K>oo,ext 


= 171.86 W/K-m, 


Koc,MC = 


172.17 ± 7.07 W/K • 


m 
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TABLE V: Thermal conductivity m, standard deviation of thermal resistivity a%, and output 
temperature Tj obtained at different sample length Lj (i = 1, 2, 8) but fixed input temperature 
of 500 K, heat flux of 0.0010 (eV/ps-A 2 ), and averaging time of 40.0 ns without thermal expansions. 



i 


1 


2 


3 


4 


5 6 


7 


8 


Li(cells) 


150 


200 


250 


300 


350 400 


450 


500 


Ki(W/K-m) 


28.72 


34.93 


39.80 


44.89 


50.82 54.79 


57.91 


62.86 


(Ti{K • m/W) 


0.000675 


0.000521 


0.000392 


0.000417 


0.000308 0.000286 


0.000261 


0.000243 


Ti 


501.04 


501.11 


500.98 


501.17 


501.23 501.22 


501.24 


501.07 




K>oo,ext 


= 120.59 W/K-m, 


Koc,MC = 


120.88 ± 5.91 W/K • 


m 





2.0 W/K • m at 800K. As will be discussed below, the uncertainty of extrapolated thermal 
conductivity is extremely sensitive to the standard deviation of thermal resistivity obtained 
from MD simulations. 

In Fig. [ill a) we show the MD results and the fit to Eq. [2] at each temperature. While 
the length dependence predicted by Eq. [2] appears reasonable, there are some systematic 
deviations that can be seen at each temperature. In particular, the data points obtained at 
longer sample length define a steeper slope than the data points at shorter sample length. 
The effect is most notable at the highest temperature of 800 K, and least significant at 
the lowest temperature of 300 K. This deviation from linearity has not been recognized 
in previous work, but is clearly important to the theoretical study of thermal transport 
properties. Since this effect is relatively insignificant at 300 K and 500 K, the overall linear 
dependence of 1/k on 1/L seems to apply very well. The values of obtained from the 



linear fits in Fig. [TT]are 184.67 W/K ■ m at 300 K and 101.64 W/K ■ m at 500 K. We include 
the extrapolated values of K m in Tables [I| - |V| as Koo,ext- We found that the values of «:„ 
obtained using Eq. [2] were different when thermal expansion was not included. From Tables 
[TJ [IT], [TV] and |VJ it can be seen that when thermal expansion was not included, the value 
was decreased to 171.86 W/K at 300 K and increased to 120.59 W/K • m at 500 K. The 
difference became more significant at 500 K than at 300 K, consistent with a larger thermal 
expansion at higher temperatures. 

The more significant deviation of the 800 K data from the predictions of Eq. [2] merits 
further careful examination because it is not clear if Eq. [2] can be used with a good confidence 
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( a) all temperature data 



(b) expanded 800 K temperature data 




0.002 0.003 0.004 0.005 0.006 0.007 0.002 0.003 0.004 0.005 0.006 0.007 



1/sample length (1/cell) 1/sample length (1/cell) 



FIG. 11: MD results for the computed resistivity 1/Ki for different inverse sample lengths 1/Lj. 
The lines represent fits of Eq. [2] to the MD data. 

to obtain k^. One possible explanation for the nonlinear behavior is that the individual 
data points «j might not correspond to the linear-transport regime. This point was explored 



above by varying the heat flux as shown in Fig. [TO] which showed no dependence on the heat 
current J. However, that investigation was done at 300 K where the nonlinear transport is 
the least significant. Here we perform several additional calculations at 800 K, some with 
very low heat flux J. These results are shown in Fig. [TT|(b). While the nonlinear behavior 
is still present for the data obtained with a reduced heat current J, there is some obvious 
dependence of Ki on the magnitude of the current J for intermediate cell sizes (e.g. for n\ 
in the neighborhood of 300). In particular, the linear relationship for the short sample end 
(Li below 300 cells) is extended to longer sample length by the reduction of the heat current 
J. In practices, the approach to use smaller heat current is limited because it increases 
the thermal fluctuation and the associated errors. Under the assumption that the linear 
behavior as seen at smaller sample sizes could be extended to larger system sizes by further 
decreasing the heat current J (requiring significantly-increased simulation time to contain 
the error), we fit Eq. [2] with only the seven data points that correspond to small system sizes 
in Fig. [TT|(b). We obtain a value of 71.88 W/K • m. While this value of is reasonable, 
it is clear that the nonlinear behavior is an issue that calls into question the use of Eq. [2] at 
least in some cases. 

The fit to Eq|2]also includes a prediction of the slope a. The value of a can be estimated 
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as a = 8/ks ■ n ■ v, where fee is Boltzmann's constant, n is the number density, and v is the 
average acoustic phonon velocity. For GaN, the number density is about 0.087 A~ 3 , and the 
average velocity v ~ 5 x 10 3 (m/s). This leads to an estimated slope of a = 1.33 xlO -9 



m 2 • K/W. For the results obtained here and shown in Fig. 11, we find slopes of a = 1.48 
xl0~ 9 m 2 ■ K/W at T = 300K, a = 1.96 xl0~ 9 m 2 ■ K/W at T = 500K, and a = 2.79 
xlO -9 m 2 ■ K/W at T = 800 K. The 300 K and 500 K data are in reasonable agreement 
with the theoretical prediction, and the 800 K data is substantially larger. This is expected 
as 800 K is above the Debye temperature. 

Finally, we compare our results with previous atomic-scale simulations and experiments. 
Wang et al[36] used Evans' homogeneous field method to compute the thermal conductivity 
of GaN using the same SW potential. For 300 K, they obtain a value of 215 W/K ■ m, 
somewhat higher than our result of 184.67 W/K-m. In another study, Kawamura and 
coworkers]?] used Green-Kubo approach to find a value in the range 310 - 380 W/K ■ m at 
300 K, which is significantly higher than those reported here and in Ref. Their value, 
however, appears to have substantial statistical error. On the other hand, the experimental 
room-temperature thermal conductivity was reported to be 170 - 180 W/K ■ m by Asnin et 
al07J, ~ 155 W/K ■ m by Luo et al[48j, 186 - 210 W/K ■ m by Florescu et al[49J, ~ 250 
W/K ■ m by Slack et al|50j, and ~ 220 W/K ■ m by Jezowski et al[5T]. Our predicted result 
is in the range of the experimental values. However, it is important to take note of the fact 
that the MD simulation results are purely classical, while 300 K is substantially below the 
Debye temperature of GaN. 



IV. ANALYSIS AND DISCUSSION 



A. Monte Carlo analysis of error of extrapolated thermal conductivity 

As shown in Tables |T] - [V] the results of MD simulations are a series of thermal conduc- 
tivities K%, K2, and K n and standard deviations of the thermal resistivity a\, <72, and 
a n , obtained at n different sample lengths Li, L 2 , and L n . The standard deviation of 
the thermal conductivity, e±, e n , can be estimated from the standard deviation of 

thermal resistivity using the approximate relation Eq. [6j However, the value of /too obtained 
from Eq. [2] might have a standard deviation that is substantially different from that of the 
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individual data points. Here we develop and present a Monte-Carlo analysis of the error in 
/too and its dependence on the standard deviation of thermal resistivity of each data point 
a t . 

We can assume that the thermal resistivity data obtained from a new set of MD simula- 
tions can be expressed as 

(I/O* = l/m + A (I/O (7) 

where can be viewed as the best estimate of the thermal resistivity and A (I/O is the 
deviation of the thermal resistivity for a hypothetical calculation. Assuming that A (I/O 
satisfies a Gaussian probability distribution, 

2" 



P [A (!/««)] = 7= — o ea; P 

2,TT(jf 



1 /A (I/O 



2 V a, 



(8) 



where cij is the standard deviation of resistivity. Because l/«, and <jj are known, we can 
randomly sample the (I/O* values that would be from hypothetical MD calculations. To 
do this, a random number r between zero and one is created. The deviation A (I/O that 
satisfies the probability distribution function of Eq. [8] can be obtained from 

-A(l/«i) 

p (x) dx = r (9) 

and (I/O* is then obtained from Eq. [Tj Once a set of (I/O* = 1> 2, n) values 
are sampled, a linear regression can be carried out using Eq. [2j and a sampled value is 
obtained. After generating many sets of data via this Monte-Carlo method, the standard 
deviation for be determined. 

To facilitate the examination of main factors controlling the error of the extrapolated 
thermal conductivity, we assume that the best estimate of the thermal conductivity (O at 
sample length (O can be described by Eq. [2] with parameters l//Coo = 0.0054234 (K • m/W) 
and a = 2.85035 (K ■ m ■ cell/W) (note that the unit of Lj = n\ used here is number of 
cells), and the standard deviation of resistivity is constant o\ = o~2 = ••• = 0n = 0.0008 
(K • m/W). These l/re^ and a values are fairly close to the best fit of the real 300 K MD 
data shown in Table |TJ A total of 100000 sets of MC data were generated. The resulting 
values of K m were used to obtain an average Koo,mc an d its standard deviation e^MC- 

The effects of several parameters are explored, including number of samples n, the min- 
imum and the maximum sample length L\ and L n (assuming that Li, L2, and L n are 
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( ft) effect of sample numbers 



(b) effect of minimum sample length 



5 




5 10 15 20 25 

number of different samples 



(c) effect of maximum sample length 



minimum sample length (lattice cell) 
(d) effect of unextrapolated standard deviation 




300 350 400 450 500 550 

maximum sample length (lattice cell) 



0.0000 0.0003 0.0006 0.0009 0.0012 0.0015 

MD standard deviation of resistivity a; (K.-m/W) 



FIG. 12: Factors determining the standard deviation of thermal conductivity at extrapolated 
infinite sample length. 



in order), and the MD (unextrapolated) standard deviation of resistivity <7j. Figs. |l2[a) - 
l~2^d) show the extrapolated standard deviation of conductivity as a function of each of the 
parameters with all the other parameters kept constant: Li = 150, L n = 500, n = 8, and a 
= Oi = 0.0008 (K • m/W) (i = 1, 2, n). 



Fig. 12 a) shows the dependence on the number of samples n with length uniformly dis- 
tributed between L\ = 150 and L n = 500. It can be seen that e^MC decreases with the 
number of samples n, consistent with the intuition that more data improves the accuracy. 
One key observation is that the number of samples cannot be too small or the error may 



increase abruptly to very large values. Figs. [I2[b) and 12 c) show that increasing the min 



imum sample length while keeping the maximum sample length fixed, or decreasing the 
maximum sample length while keeping the minimum sample length fixed, causes an increase 
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in the error e^MC- This is expected because a narrow sample length range is associated 
with a large uncertainty in the fit. The value of e^^c can increase abruptly when the 
minimum sample length is above or the maximum sample length is below some threshold. 



Finally, Fig. 12 'd) shows that reducing the standard deviation of thermal resistivity Oi for 
MD simulations has a significant effect on the standard deviation e^MC of the extrapo- 
lated conductivity. In general, e^MC is nearly proportional to Oi in the range Ui < 0.0012 
K ■ m/W. However, there exists a threshold <7j value between 0.0013 K ■ m/W and 0.0014 
K • m/W, above which e^^tc abruptly changes to very large values. The abrupt increases 



of the error seen in Figs. 12 a) - [121(d) occur because occasionally zero or unphysical (nega- 



tive) 1 / values are sampled. The finding revealed in Fig. [I2](d) indicates that in order to 
control error of extrapolated conductivity to below 10 W/K • m, the MD error of resistivity 
should be less than 0.0003 K ■ m/W (corresponding to < 1.0 W/K ■ m error in conductivity). 
Often, higher accuracy is required if the number of different sample length or the difference 
between minimum and maximum sample length are not as large as analyzed here. To our 
knowledge, the errors of most published MD data do not satisfy this requirement. As a 
result, significant differences exist among the published data. 

To facilitate the parameter study, the analysis above assumes that the best estimate of 
the thermal conductivities Hi satisfies a fitted curve and the standard deviation of resistivity 
<7j equals a constant. As discussed above, real and <7j data have been obtained from the 
MD simulations, Table |T]-|Vj These data can be directly used in the Monte Carlo approach 
to estimate the MC averaged Koo,mc an d the standard deviation e^^ic- The best values 
(i.e. using the largest possible averaging time) obtained for Koo,mc and e^Afc are shown in 
Table [I] - |Vj Figs. 13 - 15 show Koo,a/c and e^MC as a function of averaging time for the 
300 K, 500 K, and 800 K simulations. The 800 K results use only the 7 MD simulations 



in the short sample length range, Fig. [IT[b). It can be seen from Figs. 13 and 14 that 
the conductivity reached a converged value with a relatively small standard deviation (< 
10 W/K • m) at 300 K and 500 K provided that the averaging time is 10 ns or above. In 
sharp contrast, the 800 K data shown in Figs. 15 demonstrate a failure to obtain Koo,mc 
for 24 ns or less. This failure results from the occasional samples of zero or negative 1 / 
values during the MC analysis, leading to divergent behavior for Koo.a/c*- Converged values 
of ^oo,mc with relatively small deviation (e.g. between 10 W/K • m and 15 W/K ■ m) can 
be achieved for an averaging time exceeding 35 ns. The significantly increased difficulties 
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(a) MC-extrapolated thermal conductivity 
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(b) MC-derived standard deviation 
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FIG. 13: MC-averaged thermal conductivity Koo,MC and its standard deviation €oo,MC as a function 
of averaging time at 300 K. 



(a) MC-extrapolated thermal conductivity 
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(b) MC-derived standard deviation 
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FIG. 14: MC-averaged thermal conductivity Kqo, mc an d its standard deviation e^^MC as a function 
of averaging time at 500 K. 



in extrapolating the 800 K data are a consequence of both reduced sample length range, 
Fig. |l~2^c), and increased thermal fluctuations at higher temperature, Fig. [l2^d) . 

To explore the origin of the divergence in the average value Koo, mc found in the 800 K MD 
results, we show in Fig. fl6Fa) the probability distributions for obtained for averaging 



times of 8 ns and 40 ns. It can be seen from Fig. 16 ^a) that for 8 ns averaging time at 



800K, the distribution deviates strongly from a Gaussian and exhibits a long tail to very 
large values of k^. It is exactly this long tail that leads to diverging values of Koo,mc an d 
very large standard deviations e^MC- Increasing the averaging time to 40 ns results in a 
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(a) MC-extrapolated thermal conductivity (b) MC -derived standard deviation 
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FIG. 15: MC-averaged thermal conductivity Koo,mc an d its standard deviation e^^C as a function 
of averaging time at 800 K. 



(a) 800 K data (b) 500 K data 
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FIG. 16: Probability density distribution for Koo- 

distribution that is more nearly Gaussian and does not have the tail for extremely large 
values. While the distribution is still fairly wide for 40 ns of averaging, it is apparent that a 
value within 10 W/m ■ K of the best (peak) result for k m will occur with a high probability. 
On the other hand, it is noticeable that for the short averaging time of 8 ns, the peak in 
the probability density is significantly displaced from that obtained from the long averaging 
time of 40 ns. For comparison, a distribution obtained from the 500 K MD simulation at 
an averaging time of 40 ns is shown in Fig. [T6[b). The distribution at 500 K is an almost 
perfectly Gaussian and sharply peaked at the mean k CO] mc value. As a result, the long time 
simulation at 500 K results in a highly accurate estimate of Koo,mc- 
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While the averaged obtained at short averaging times can give unphysical results, it 
is still apparent that the probability density near the best estimate of at the same 

time be quite large. This suggests that there is still a significant probability of obtaining 
a quite reasonable value of from a single set of simulations even when the average 
of the distribution is unphysical. The problem emerges from the small but not insignificant 
possibility of obtaining an extremely large or even divergent value of k^. This might explain 
why many published simulation results using the direct method and Eq. [2] to determine n^, 
obtain reasonable results even with much shorter simulation times. What we have uncovered 
in the present study is that without adequate averaging time and an adequate data set 
including many different system lengths Li there is a finite chance of obtaining a result for 
/too that deviates significantly from the correct value and may be extremely large. 

To further demonstrate the origin of the long tail shown in Fig. fl6~|(a), we calculated 
the probability distribution for the MC extrapolated 1/kqo values. The results at 800 K 



for averaging times of 8 ns and 40 ns are shown in Fig. 17 As expected, the distribution 
for the long averaging time is relatively sharply peaked and the significant 1 / values are 
larger than zero. For the 8 ns averaging time, the distribution is still Gaussian but much 
broader. More importantly, the distribution extends to low side of the peak point, and 
even reaches zero and the negative range. The small 1 / correspond to very large values 



for /too, resulting in the tail seen in Fig. 16 a). In computing k^mg as the average of the 
distribution, these very large values of lead to the failure to obtain a converged value in 
some instances. 



B. Simulation time to reach steady state 

The results presented above show that to be assured of accurate results, extremely long 
averaging times are required. One possible approach to decrease the overall simulation time 
is to initialize the temperature distribution to be close to the expected steady-state distri- 
bution. To explore this, simulations were carried out for both a uniform and a linear initial 
temperature distribution at an average sample temperature of 300 K, a sample length of 500 
cells, a heat flux of 0.0015 eV/ps ■ A 2 , and a source width of 60 A. The temperature data 



averaged between 0.95 ns and 1.0 ns are shown in Fig. [I8[a) and Fig. 18 b) respectively for 



the uniform and the linear initial temperature distribution. For comparison, the initial tem- 
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FIG. 17: Probability density distribution of 1/kqq, k p is the peak value of the k distribution. 



perature profiles are shown using the black lines, and the steady-state temperature profiles 
(taken as the one that is fully averaged between 4.0 ns and 24.0 ns, see Fig. [6]) are shown 
using the gray lines. It can be seen that at an absolute time of 1.0 ns, the temperature 
profiles averaged between 0.95 ns and 1.0 ns for the uniform and linear initial temperature 
distributions are very similar, and the memory of initial temperature distribution has been 
lost. This means that the acceleration of the simulation by better initial temperature distri- 
bution does not play a major role in the computation cost. The time required to achieve the 
steady-state during heat diffusion can be estimated from diffusion distance d and diffusivity 
D using 



D 



(10) 



where d = L/2, D = k/C v , and C v is volumetric heat capacity. In classical MD simulations, 
heat capacity C v = 3p-ks, where p is material density. For the sample simulation presented 



in Fig. 18, « = 90.76 W/K • m, d = 1302 A, p = 0.0868 atoms/A 3 , resulting in r w 0.67 ns. 
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(a) uniform initial temperature distribution 
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(b) linear initial temperature distribution 
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FIG. 18: Evolution of temperature profiles from different initial temperature distributions. 



This is certainly consistent with the results shown in Fig. 18 



C. Green-Kubo thermal conductivity calculations 



To corroborate the results of the direct method, the Green-Kubo method was used to 
compute the thermal conductivity at 500 K and 800 K. In the Green-Kubo method, the 
thermal conductivity is usually expressed as[ 

1 



lim 

T— »0O 



{Ji{t)Jj{G))dt 



(11) 



Jo 

where Q is the system volume and Ji(t) is the ith component of the thermal current. Nu- 
merically, this integral of the autocorrelation can be estimated by 



KIT 



MAt) 



At Jj(m + n)Jj(n) 



nk B T 2 



N — m 



M < N 



(12) 



m=l n=l 

The calculations were carried out at a time step size of At = 0.32 fs. It is known that 
the Green-Kubo method converges slowly especially at high temperatures. As a result, four 
separate simulations, each with 10 7 time steps, were performed at 500 K, and 24 separate 
simulations, each with 9 x 10 6 time steps, were performed at 800 K. These represent about 
13 ns and 69 ns total simulation time at the two temperatures, thereby effectively reducing 
the thermal fluctuations. 



We used the bootstrap algorithm [52] to compute the thermal conductivity from Eq. 12 
The bootstrap procedure involves choosing random samples from a given data with replace- 
ment in order to estimate statistics. In our case, we used the method to calculate confidence 
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intervals for the integral of the heat flux autocorrelation Eq. 11 corresponding to the outer 



sum on m in Eq. 12 by using random samples of the inner sum on n in Eq. 12 representing 
the mean autocorrelation at index m. Since the ordinary bootstrap assumes that data are 
independent and our data are strongly correlated, we used a variant of the bootstrap that 
employs block averaging. We have about 3200 ps of data for each simulation. We have 
created 100 blocks of about 32 ps from these data by averaging them. Then we applied 
the bootstrap on these blocks [52| p. 396] with the assumption that correlation between 
observations is strongest within a block and relatively weak between blocks. 

To illustrate the results, the averaged current-current correlation functions along the 



[1100] direction are shown in Fig. 19 for both 500 K and 800 K temperatures. The results 



obtained along the other two directions are similar. It is interesting to see from Fig. 19 



that the current-current correlation function exhibits two fast oscillations, corresponding 
well with the two optical modes. The amplitude decays rapidly as the relaxation time is 
increased. Specifically, the current-current correlation function has reached near-zero when 
the relaxation time is increased to 30 ps. 

The thermal conductivities in all the three coordinate directions [1120], [1 100] and [0001] 



are shown in Fig. 20 and 21 respectively for the 500 K and 800 K temperatures. For 500 



K, Fig. 20 shows that the approximately converged thermal conductivity at a relaxation 
time of 50-60 ps is around 77 ± 14 W/K ■ m in the [1100] direction, 88 ± 16 W/K • m in 



the [1120] direction, and 98 ± 16 W/K • m in the [0001] direction. For 800 K, Fig. [21] shows 
that the approximately converged thermal conductivity is around 58 ± 9 W/K ■ m in the 
[1100] direction, 58 ± 8 W/K • m in the [1120] direction, and 71 ± 9 W/K ■ m in the [0001] 
direction. 

Despite the large scale of the calculations especially for the 800 K case, the results are still 
associated with relatively significant statistical errors. This is not unusual for the Green- 
Kubo method. Nonetheless, the results obtained from the Green-Kubo method corroborate 
well with those from the direct method. 



V. CONCLUSIONS 

Large-scale parallel computer simulations have been carried out to identify sources of 
errors of the commonly used direct molecular dynamics method for thermal conductivity 
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time (ps) 

FIG. 19: Current-current correlation function < J(t) ■ J(0) > in the [1100] direction as a function 
of relaxation time r at 500 K and 800 K. 

calculation and ways that can significantly reduce the errors. As a case study, we investi- 
gated the thermal transport along the [0001] direction of a model GaN bulk crystal using 
a Stillinger- Weber potential[3U [35]. A Monte Carlo method has been developed to quan- 
tify the uncertainty of the thermal conductivity values that are extrapolated to an infinite 
sample length from simulated values at finite sample length. 

From simulations of extremely long duration, we have obtained results with a level of 
accuracy not previously attained. The high accuracy and careful statistical analysis of our 
results have enabled a detailed investigation of the validity of direct method calculations 
of thermal conductivity including the procedure to determine Koo from the predicted size 
dependence of Eq. [2j We have reached several important conclusions: 

(a) For the direct method, the averaging time required to produce accurate results depends 
on material, system dimension, and temperature. Smaller system cross-sectional area 
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FIG. 20: Thermal conductivity along the three coordinate directions as a function of relaxation 
time at 500 K. 

and higher temperature generally demand longer averaging time. Most of our calcula- 
tions used an averaging time of 40 ns. This is significantly longer than that commonly 
used in the literature; consequently, our results are considerably more accurate. 

(b) With tightly controlled error relative to existing work, our results provide strong ev- 
idence that within an appropriate parameter range, model parameters such as cross- 
sectional area, heat source width, and heat flux do not affect the results of thermal 
conductivity. 

(c) Computed thermal conductivities «j for different system lengths Lj show in most cases 
good agreement with the length dependence predicted by Eq. [2j However, clear devi- 
ations from Eq. [2] are observed especially at high temperatures, large sample lengths, 
and relatively large thermal currents J. 
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FIG. 21: Thermal conductivity along the three coordinate directions as a function of relaxation 
time at 800 K. 

(d) Deviations between the predictions of Eq. [2] and MD results appear to be due to 
nonlinear transport effects. Reducing the heat current J leads to linear behavior for 
longer system sizes. 

(e) Monte Carlo analysis shows that the extrapolated conductivity is extremely sensi- 
tive to the quality of the data set, including the standard deviation of each simulation, 
number of system lengths, and the minimum and maximum sample lengths. Occasion- 
ally extremely large or unphysical values of be obtained when an insufficient 
data set is used to fit Eq. [2} 

(f) The analysis performed here suggests that very long simulations ~ 20 - 40 ns might 
be required to assure accurate results. Probability distributions from MC results, 
however, indicate that previous studies based on much shorter simulation times ~ 1 - 
2 ns have a significant likelihood of producing reasonable results. This seems to explain 
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why Eq. [2] has occasionally failed to give reasonable results, but yet often yields values 
for in good agreement with other simulation approaches and experiment. 

We predict that the [0001] thermal conductivity of GaN bulk crystal is 185 W/K • m at 
300 K, 102 W/K • m at 500 K, and 74 W/K • m at 800 K. These compare really well with 
98 W/K • m at 500 K and 71 W/K • m at 800 K predicted from the Green-Kubo method. 

We believe that the observed nonlinear behavior in the dependence of 1/k on 1/L is 
perhaps the most significant outcome of this work. However, we do not necessarily con- 
clude that Eq. [2] is not able to correctly predict the length dependence. In particular, the 
occurrence of the nonlinear behavior seems to depend on the magnitude of the thermal cur- 
rent J, which suggests that it should be possible, if not always practical, to compute for 
longer system lengths Lj with a smaller current J and recover the linear behavior. However, 
we have also found that uncovering the nonlinear dependence requires extremely accurate 
calculations which can only result from averaging over rather long simulation times. 

In summary, the results presented here suggest that, in some circumstances, considerably 
more care might be required to determine 1/k,oc using the direct method than has often been 
exercised. Even fairly small standard deviations in each calculated Ki can lead to rather large 
errors when using Eq. [2] to determine Kqo. To overcome these difficulties, extremely long 
simulation times (> 20 ns) might be required. Perhaps more importantly, very accurate 
data might be required to identify nonlinearities which can lead to unphysical values of 
when fitting data using Eq. [2j 
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